function gout=freqfilter(winname,bw,varargin)
%FREQFILTER Construct filter in frequency domain
%   Usage:   g=freqfilter(winname,bw);
%            g=freqfilter(winname,bw,fc);
%
%   Input parameters:
%      winname  : Name of prototype
%      bw       : Effective support length of the prototype
%      fc       : Center frequency
%   
%   `freqfilter(winname,bw)` creates a full-length frequency response
%   filter. The parameter *winname* specifies the shape of the frequency
%   response. For accepted shape please see |freqwin|. *bw* defines a
%   -6dB bandwidth of the filter in normalized frequencies.
%   
%   `freqfilter(winname,bw,fc)` constructs a filter with a centre
%   frequency of *fc* measured in normalized frequencies.
%
%   If one of the inputs is a vector, the output will be a cell array
%   with one entry in the cell array for each element in the vector. If
%   more input are vectors, they must have the same size and shape and the
%   the filters will be generated by stepping through the vectors. This
%   is a quick way to create filters for |filterbank| and |ufilterbank|.
%   
%   `freqfilter` accepts the following optional parameters:
%
%     'fs',fs     If the sampling frequency *fs* is specified then the 
%                 bandwidth *bw* and the centre frequency *fc* are 
%                 specified in Hz.
%
%     'complex'   Make the filter complex valued if the centre frequency
%                 is non-zero. This is the default.
%
%     'real'      Make the filter real-valued if the centre frequency
%                 is non-zero.
%
%     'delay',d   Set the delay of the filter. Default value is zero.
%
%     'scal',s    Scale the filter by the constant *s*. This can be
%                 useful to equalize channels in a filter bank.
%
%     'pedantic'  Force window frequency offset (g.foff) to a subsample 
%                 precision by a subsample shift of the window.
%

% Authors: Nicki Holighaus & Zdenek Prusa
% Date: September 15, 2016 

if ~iscell(winname), winname = {winname}; end

% Define initial value for flags and key/value pairs.
definput.import={'normalize'};
definput.importdefaults={'energy'};
definput.keyvals.delay=0;
definput.keyvals.fc=0;
definput.keyvals.fs=2;
%definput.keyvals.order=4;
definput.keyvals.scal=1;
definput.keyvals.min_win=1;
%definput.keyvals.trunc_at=10^(-5);
definput.keyvals.bwtruncmul = 4;
definput.flags.pedantic = {'pedantic','nopedantic'};
definput.flags.real={'complex','real'};

[flags,kv]=ltfatarghelper({'fc'},definput,varargin);

[bw,kv.fc,kv.delay,kv.scal]=scalardistribute(bw,kv.fc,kv.delay,kv.scal);

% Sanitize
kv.fc=modcent(2*kv.fc/kv.fs,2);

Lw = @(L,bw) min(ceil(bw*kv.bwtruncmul*L/kv.fs),L);
    
fsRestricted = @(L,bw) kv.fs/L*Lw(L,bw);
if flags.pedantic
    fc_offset = @(L,fc) L/2*fc-round(L/2*fc);
else
    fc_offset = @(L,fc) 0;
end

Nfilt = numel(bw);
gout = cell(Nfilt,1);
for ii=1:Nfilt
    g=struct();
    g.foff=@(L) round(L/2*kv.fc(ii)) - floor(Lw(L,bw(ii))/2); 
    
    if flags.do_1 || flags.do_area 
        g.H=@(L) fftshift(...
                 freqwin(winname,Lw(L,bw(ii)),bw(ii),fsRestricted(L,bw(ii)),...
                 'shift',fc_offset(L,kv.fc(ii)))...
                 )*kv.scal(ii)*L;                
    end
    
    if  flags.do_2 || flags.do_energy
        g.H=@(L) fftshift(...
                 freqwin(winname,Lw(L,bw(ii)),bw(ii),fsRestricted(L,bw(ii)),...
                 'shift',fc_offset(L,kv.fc(ii)))...  
                 )*kv.scal(ii)*sqrt(L);                        
    end
        
    if flags.do_inf || flags.do_peak
        g.H=@(L) fftshift(...
                 freqwin(winname,Lw(L,bw(ii)),bw(ii),fsRestricted(L,bw(ii)),...
                 'shift',fc_offset(L,kv.fc(ii)))... 
                 )*kv.scal(ii);           
    end
          
    g.realonly=flags.do_real;
    g.delay=kv.delay(ii);
    g.fs=kv.fs;
    gout{ii}=g;
end

if Nfilt==1
    gout=g;
end;
